REDPIC

Tutotial

V. Fedorov

In [1]:
import redpic as rp
import kenv as kv
import numpy as np
import holoviews as hv
hv.extension('matplotlib')

import warnings
warnings.filterwarnings('ignore')
In [2]:
rp.__version__
Out[2]:
'1.0.0.0'
In [3]:
kv.__version__
Out[3]:
'0.2.2.0'

Some plotting options:

In [4]:
%output size=150 backend='matplotlib' fig='png' dpi=200
%opts Curve Scatter [aspect=3 show_grid=True]
%opts Curve (linewidth=1 alpha=0.7 color='blue')
%opts Scatter (alpha=0.7 s=0.5)

Define accelerator beamline parameters:

In [5]:
acc = rp.Accelerator(0.7, 7.7, 0.01)
In [6]:
#              Unique name,  z-position [m],  Ez [MV/m],  Ez(z) profile
acc.add_accel('Acc. 1',      4.096,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 2',      5.944,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 3',      6.796,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 4',      8.644,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 5',      9.496,          -1.1,         'Ez.dat')
In [7]:
#                 Unique name,  z-position [m],  Bz [T],  Bz(z) profile
acc.add_solenoid('Sol. 1',      0.45,           -0.058,   'Bz.dat')
acc.add_solenoid('Sol. 2',      0.957,           0.039,   'Bz.dat')
acc.add_solenoid('Sol. 3',      2.107,           0.025,   'Bz.dat')
acc.add_solenoid('Sol. 4',      2.907,           0.044,   'Bz.dat')
acc.add_solenoid('Sol. 5',      3.670,           0.04,    'Bz.dat')
acc.add_solenoid('Sol. 6',      4.570,           0.0595,  'Bz.dat')
acc.add_solenoid('Sol. 7',      5.470,           0.059,   'Bz.dat')
acc.add_solenoid('Sol. 8',      6.370,           0.060,   'Bz.dat')
acc.add_solenoid('Sol. 9',      7.270,           0.065,   'Bz.dat')
acc.add_solenoid('Sol. 10',     8.170,           0.065,   'Bz.dat')
acc.add_solenoid('Sol. 11',     9.070,           0.0655,  'Bz.dat')
acc.add_solenoid('Sol. 12',     9.970,           0.075,   'Bz.dat')
In [8]:
acc.compile()
In [9]:
dim_z  = hv.Dimension('z',  unit='m')
dim_Ez = hv.Dimension('Ez', unit='MV/m', label='$E_z$')
dim_Bz = hv.Dimension('Bz', unit='Gs', label='$B_z$')
In [10]:
z  = acc.z
z_Ez = hv.Curve((z, acc.Ez(z)), kdims=dim_z, vdims=dim_Ez)
z_Bz = hv.Curve((z, acc.Bz(z)*1e4), kdims=dim_z, vdims=dim_Bz)
In [11]:
(z_Ez + z_Bz).cols(1)
Out[11]:
In [12]:
print(acc)
Accelerator structure.
	Solenoids:
	[ 0.45000 m, -0.05800 T, 'Bz.dat', 'Sol. 1', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 0.95700 m, 0.03900 T, 'Bz.dat', 'Sol. 2', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 2.10700 m, 0.02500 T, 'Bz.dat', 'Sol. 3', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 2.90700 m, 0.04400 T, 'Bz.dat', 'Sol. 4', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 3.67000 m, 0.04000 T, 'Bz.dat', 'Sol. 5', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 4.57000 m, 0.05950 T, 'Bz.dat', 'Sol. 6', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 5.47000 m, 0.05900 T, 'Bz.dat', 'Sol. 7', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 6.37000 m, 0.06000 T, 'Bz.dat', 'Sol. 8', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 7.27000 m, 0.06500 T, 'Bz.dat', 'Sol. 9', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 8.17000 m, 0.06500 T, 'Bz.dat', 'Sol. 10', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 9.07000 m, 0.06550 T, 'Bz.dat', 'Sol. 11', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 9.97000 m, 0.07500 T, 'Bz.dat', 'Sol. 12', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	Accelerating modules:
	[ 4.09600 m, -1.10000 T, 'Ez.dat', 'Acc. 1', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 5.94400 m, -1.10000 T, 'Ez.dat', 'Acc. 2', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 6.79600 m, -1.10000 T, 'Ez.dat', 'Acc. 3', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 8.64400 m, -1.10000 T, 'Ez.dat', 'Acc. 4', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	[ 9.49600 m, -1.10000 T, 'Ez.dat', 'Acc. 5', 0.00000 m, 0.00000 rad, 0.00000 m, 0.00000 rad],
	Quadrupoles:
	Correctors x:
	Correctors y:

Define beam parameters:

In [13]:
beam = rp.Beam(
    type=rp.electron, 
    energy = 1.32,          # MeV
    current = 0.5e3,  # A
    radius_x = 48e-3, # initial r (m)
    radius_y = 48e-3, # initial r (m)
    radius_z = 3.5,
    radius_xp = 2*35.0e-3,     # initial r' (rad)
    radius_yp = 2*35.0e-3,     # initial r' (rad)
    x  = 0.0e-3,   # horizontal centroid position (m)
    xp = 0.0e-3,     # horizontal centroid angle (rad)
    y = 0,          # vertical centroid position (m)
    normalized_emittance = 200e-6*np.pi*3.58) # m*rad
In [14]:
beam.generate('KV', 80_000)
In [15]:
beam.df
Out[15]:
x y z px py pz
0 -0.025228 -0.032257 0.225161 -0.068319 -0.091167 1.771276
1 0.036267 0.029977 -3.202804 0.105126 0.085247 1.743588
2 -0.010887 -0.030046 -0.138010 -0.037309 -0.094209 1.760767
3 0.045689 0.007946 1.965987 0.134787 0.020621 1.763194
4 -0.036430 -0.020652 -1.198027 -0.100615 -0.060292 1.755270
... ... ... ... ... ... ...
79995 -0.007845 -0.020223 2.253929 -0.015449 -0.051869 1.750142
79996 0.031381 0.002904 -3.167089 0.098358 0.003079 1.746981
79997 0.025505 -0.035000 0.278731 0.079251 -0.103145 1.767815
79998 -0.006563 0.027261 3.025015 -0.018879 0.070151 1.755272
79999 0.025116 -0.018896 0.643528 0.069158 -0.047418 1.744894

80000 rows × 6 columns

In [16]:
dim_x = hv.Dimension('x', unit='m', range=(-0.1, 0.1))
dim_y = hv.Dimension('y', unit='m', range=(-0.1, 0.1))
dim_z = hv.Dimension('z', unit='m', range=(acc.z_start, acc.z_stop))
dim_px = hv.Dimension('px', unit='MeV/c', label='$p_x$')
dim_py = hv.Dimension('py', unit='MeV/c', label='$p_y$')
In [17]:
beam_x_y = hv.Scatter(beam.df, kdims=[dim_x, dim_y])
beam_z_x = hv.Scatter(beam.df, kdims=[dim_z, dim_x])
beam_x_px = hv.Scatter(beam.df, kdims=[dim_x, dim_px])
beam_y_py = hv.Scatter(beam.df, kdims=[dim_y, dim_py])
In [18]:
(beam_x_y + beam_z_x + beam_x_px + beam_y_py).cols(2)
Out[18]:
In [19]:
print(beam)
Beam parameters:
	Type	electron
	Distribution	KV
	Particles	80000
	Current	500 A
	Energy	1.320 MeV
	Total momentum	1.758 MeV/c
	Rel. factor	3.583
	Radius x	48.0 mm
	Radius y	48.0 mm
	Radius z	3.5 m
	Radius x prime	70.0 mrad
	Radius y prime	70.0 mrad
	Horizontal centroid position	0.0 mm
	Vertical centroid position	0.0 mm
	Horizontal centroid angle	0.0 mrad
	Vertical centroid angle	0.0 mrad
	Normalized emittance x	2249.4 mm*mrad
	Normalized emittance y	2249.4 mm*mrad

Run simulation!

In [20]:
kv_sim = kv.Simulation(beam, acc)
kv_sim.track()
In [21]:
rp_sim = rp.Simulation(beam, acc)
rp_sim.track()
z = 7.68 m (99.7 %) 

Plot the simulation results:

In [22]:
def plot(i):
    df = rp_sim.result[i]
    kv_z_x = hv.Curve(((acc.z, kv_sim.envelope_x(acc.z))), kdims=[dim_z], vdims=[dim_x], label='kenv')*\
    hv.Curve(((acc.parameter,-kv_sim.envelope_x(acc.z))), kdims=[dim_z], vdims=[dim_x])
    rp_z_x = hv.Scatter(df, kdims=[dim_z, dim_x], label='redpic')
    return rp_z_x*kv_z_x
In [23]:
items = [(i, plot(i)) for i in list(rp_sim.result.keys())]

hv.HoloMap(items, kdims = ['z']).collate()
Out[23]:
In [ ]: